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^ \ Abstract: We present a brief description of the recently developed Fortran 

lyQ \ code TERAD91 for semi-analytical calculations of the double differential 

P ' cross sections of NC and CC deep inelastic electron-proton scattering and of 

some related observables. The code is mainly intended for the calculations 
K> , at HERA or other ep-coUiders but may be also used for similar processes like 

H \ muon-proton scattering in fixed-target experiments. 



1 Program Summary 

Title of the program: TERAD91 

Version: version 2.10, 3 Oct. 1991 

Computer: IBM-3090, DESY 

Program language: FORTRAN-77 

Number of program lines: about 11,000 

Other programs used: DIVDIF, interpolation routine from the CERN li- 

brary KERNLIB. 



External files needed: None. 

Method of solution: Numerical integration. 

Typical running time: Depending on the actual requirements like parame- 

ters, requested contributions, see text. 



2 Introduction 

TERAD91 is a software product by the Dubna-Zeuthen Radiative Corrections Group 
(DZRCG). It is intended for calculations of Deep Inelastic Scattering (DIS) cross sections 
and related observables (cross section ratios, asymmetries etc.) at HERA. Its typical 
output is the double differential cross section of DIS with account of QED and electroweak 
(EW) radiative corrections as chosen by flags: 

(fa 



dxdy'' 



(in nanobarn) with x and y being some scaling variables - either leptonic or mixed or 
hadronic. For a brief description of the basic theoretical issues and definitions see [1]. 

TERAD91 was created during the 1990-91 Workshop on Physics at HERA. It accu- 
mulates about 15 years experience in the field of DIS collected by the DZRCG. Formally, 
it originates from four different codes: TERAD90 [2], DISEPNC [3], DISEPCC [4] 
and DIZET [5]. 

TERAD90 was its prototype, created just before the Workshop. Its roots date back 
till the mid seventies [6, 7]. In the eighties, the TERAD approach was used for the analysis 
of muon DIS data by the BCDMS experiment [8, 9]. 

DISEPNC and DISEPCC are codes developed in the late eighties for the QPM 
calculations of DIS at HERA [3, 4]. A similar approach was used earlier for the analysis 
of neutrino DIS data by the CDHSW [10] and CHARM-I [11] experiments. 

DIZET is our electroweak library, used before for LEP-I physics [5]. It allows to 
calculate Ar, and the complete 0{a) corrections, with inclusion of some leading higher 
order terms connected with the top-quark, to Z, W widths, and to the weak form factors 
for cross sections in the s- and t-channels. These form factors depend on the kinematics 
(s, t,-u-invariants), and on the fermion type of the scattering particles. 

Some theoretical basis of the TERAD approach is presented in the contributions of 
the working group on radiative corrections in these proceedings [1]. Here we concentrate 
on the description of the code itself assuming that the reader is aware of the relevant 
physical problems and of the terminology. We present a description of its MAIN routine 
and of SUBROUTINE SET where all flags are collected which should be set by a user. 



3 MAIN program of TERAD91 

The MAIN routine of TERAD91 consists of a CALL of SUBROUTINE SET which sets various 
user flags followed by a CALL activating one of the calculational chains, which computes 
some observable (e.g. cross section) or some other physical quantity (e.g. an electroweak 
form factor). Normally only one chain is active, CALLs to others are commented. We 
recommend to work only with one chain at a time. Otherwise the running would be 
too cumbersome as some of the calculations are time comsuming and default outputs 
are lengthy. Moreover, a destructive cross influence of chains is not excluded. Below we 
briefly describe each chain. 
DXYLEP 

SUBROUTINE DXYLEP calculates double differential cross sections in terms of usual leptonic 
scaling variables x,y. For leptonic type variables, TERAD91 is mostly elaborated. It 
calculates both NC and CC cross sections in both TERAD and DISEP approaches (we 
will refer to them as two different branches). For this reason the output of DXYLEP 
contains a lot of cross sections (Born, corrected by EW and/or QED) and corrections 
calculated by both branches. The resulting cross section and correction for NC reaction 
is a mixture of numbers, calculated by both TERAD and DISEP branches. Because of 
two-fold numerical integrations, involved in the calculation of cross sections in this case, 
the running with DXYLEP is rather slow and takes about 70 minutes for a calculation in 
96 kinematical points (default option) on the DESY IBM. 
AXYLEP 

SUBROUTINE AXYLEP calculates a double differential polarization asymmetry in leptonic 
variables x, y. All calculations again are performed by both branches; since now one has 
to calculate each cross section twice, for two polarizations, the calculational time needed 
is approximately doubled. When using this chain, one should have in mind that ini- 
tialization of polarizations is done inside AXYLEP. Also electroweak parameters could be 
re-initialized there in order to investigate the sensitivity of the asymmetry to the EW 
parameters. 
DXYMIX 

SUBROUTINE DXYMIX calculates double differential cross sections in terms of mixed scaling 
variables Xm,ym- 
DXYHAD 

SUBROUTINE DXYHAD does the same in terms of hadronic variables Xh,yh (for their defini- 
tions see [1, 2]). 

For mixed and hadronic variables, the TERAD91 package contains only the TERAD 
approach for only the NC reaction. For this reason the outputs of DXYMIX and DXYHAD 
are very short: only the NC Born cross section and total correction are printed for each 
chosen {x, y) pair. As only one numerical integration is left for the calculation of the cross 
section, these chains are much faster as compared to the first two. Each of them require 
less than 1 minute for the calculation in 96 default kinematical points. 

The three subroutines DXYLEP, DXYMIX and DXYHAD are the main interfaces to the core 
of the TERAD91 code - which realizes some formulae describing differential DIS cross 
sections with account of QED and EW Radiative Corrections (EWRC) in some scaling 
variables x,y , see Section 3.2.1 of [1] and [2]. But, TERAD91 can easily also be used to 
calculate many other observables, related to differential cross sections, and the discussed 
above AXYLEP gives an example of a user interface for the calculation of another observable. 



namely polarization asymmetry A. Similar to DXYLEP, DXYMIX and DXYHAD, AXYLEP also 
returns a double differential in x and y quantity. 

The three subsequent interfaces exhibit examples of the calculation of integrated quan- 
tities: cross section ratio R = NC/CC (RMINl and RMIN2) and again the polarization 
asymmetry A (AMINl) in a kinematical region restricted by some cuts. All these sub- 
routines deal with definitions of cross sections in terms of leptonic variables and use the 
DISEP approach. In RMINl, 2 and AMINl, corresponding observables are constructed out 
of the Improved Born Approximated (IBA) DIS cross sections, i.e. with EW but without 
QED radiative corrections. This is done to study their sensitivities to the variation of the 
EW parameters. Due to additional integrations over the restricted by cuts kinematical 
region, these interfaces are rather slow. It is another reason why we use the DISEP ap- 
proach there and restrict ourselves to the IBA. 
RMINl 
AMINl 

SUBROUTINES RMINl and AMINl calculate R and A with cuts in x,y and Q^. 
RMIN2 

SUBROUTINE RMIN2 calculates R with a cut in p'l and Ql- 

All of them return R and A as functions of varying top-quark and Higgs masses. Similar 
interfaces can be readily written by a user to calculate other integrated observables. 

The three last interfaces are intended for a calculation of purely EW effects. They are 
indeed interfaces only to the EW library DIZET. They were used to produce some Ta- 
bles and Figures for the Workshop proceedings and are left in the TERAD91 as further 
examples of interfacing. 
DIZETF 

SUBROUTINE DIZETF returns a Table for the static EW parameters, Mw and sin^ 9w- 
FORMES 
FORMFT 

SUBROUTINES FORMFS and FORMFT return some tables illustrating the s- and t-dependences 
of EW form factors. We will not describe them in detail here. 



4 SUBROUTINE SET 

In this SUBROUTINE we collected all fiags used to initialize the chains. The main interfaces 
DXYLEP, DXYMIX and DXYHAD are initialized by SET completely. Other interfaces may 
contain some additional initialization inside them; SET prints also all fiags. 

In TERAD91 there are fiags which act in both the TERAD and DISEP branches; 
but there are also fiags inherent to only one of the branch. The notion of default fiags 
is not supported in TERAD91. Flag setting should be controlled by the user for each 
specific run. 



4.1 TERAD - flags 

MODEL 

M0DEL=1: quark distributions are used for the calculation of the structure functions. 



(presently the Duke and Owens [12] parametrization) , 

MODEL=0: in principle any parametrization can be used for the eight structure functions 

Fi, F2, Gi, G2, G3, Hi, H2, H3. The user can substitute instead of using SUBROUTINE SFFl 

(for Fi) etc. any parametrization available to him/her. Now we use a merge of Brasse 

[13], Stein [14] and Duke-Owens [12] fits (see [15] for more details). 

IVAR 

IVAR=0: the quark distributions are not modified in the low Q^ region, 

IVAR=1,2: they are made vanishing at low Q^ through two different modifications. They, 

as well as related physics, are described in [16]. 

ITERAD 

ITERAD=0: the quark distributions and aQED (if running, see below the description of flag 

IVPOL) in leptonic bremsstrahlung are artificially chosen to depend on Q^ as defined from 

lepton momenta. This option was implemented to have a possibility to compare results 

calculated by both TERAD and DISEP branches. 

ITERAD=1: they depend on Ql as defined from the hadronic momenta {Q\ is integration 

variable). This option must be used for realistic calculations. 

IHCUT 

IHCUT=0: no cut on the hadronic final state, 

IHCUT=1: rejects hadronic final states with Q| < TCUT (GeV^) and invariant hadronic 

mass (photon not included) W^ < AMF2CT (GeV^). In this case the two desired cut values 

TCUT and AMF2CT are set. 

IGCUT 

IGCUT=0: no cut on photon variables, 

IGCUT=1: rejects all final states with a photon energy E^ > Ecut{GeV) and with photon 

angle in the interval {9cut2,0cuti), angles in radians (0,7r) with respect to the electron 

beam. 

These two cuts (IHCUT and IGCUT) are implemented only in DXYLEP. 

IWEAK 

This flag governs how EWRC are applied in the TERAD branch. 

IWEAK=0: they are applied only via sin^ 9w, 

IWEAK=1: they are applied via the full set of the EW form factors; a QPM treatment is 

used for the structure functions Gi and Hi in this case. 

I ZERO 

a flag to swich on/off a time comsuming CALL of the TERAD branch in DXYLEP, 

IZERO=0: it is switched off; then the numbers of the TERAD branch are set to zero, 

IZER0=1: it is switched on. 



4.2 DISEP - flags 

I NORM 

INORM=0: the nine different QED corrections for the NC reaction [3] are all normalized 

by the 7-exchange Born cross section, 

IN0RM=1: they are normalized by the complete Born cross section, 7 + Z-exchange, 

IN0RM=2: each QED correction is normalized by its corresponding Born contribution, 

three with 7- exchange, three with 7Z-interference and three with Z-exchange (three 

means: lepton, lepton-quark interference, and quark bremsstrahlung). 

lEWRC 



IEWRC=0: no weak loop corrections are taken into account, 

IEWRC=1: they are taken into account. 

IBOXF 

IB0XF=1: the contributions from ZZ- and W^W^-boxes are included, 

IBOXF=0: they are not included. 

(The IBOXF flag acts only if lEWRC 7^ 0.) 

ICONV 

IC0NV=1: weak form factors are running, 

ICONV=0: they do not run. The recommended value of this flag is 0. 

(The ICONV flag is active only for IEWRC=1). 

IQ20 

IQ20=0: the light quark masses will be replaced by a Q\^ which is recommended to be 

chosen equal to the Q\ which enters the QCD evolution equations, 

IQ20=1: then the light quark masses are used. 

The first choice is recommended, the second is supported for the sake of comparison. 



4.3 Flags common both to TERAD and DISEP branches 

I EXP 

IEXP=1: the soft photon part of leptonic QED corrections is exponentiated, 

IEXP=0: it is not exponentiated. 

IVPOL 

IVPOL=0: no running of agsD, olqed = a, 

IVP0L=1: running agEv- 

IMOMS 

This flag defines the calculational scheme used, i.e. which EW parameters are considered 

as input quantities and which one is calculated (iterated). 

IM0MS=1 

IM0MS=2 

IM0MS=3 



input - a, Gfj,, Mz, rrit, Mh] Mw is iterated, 
input - a, Gfj,, Mw, rrit, Mh] Mz is iterated, 
input - a, Mw, Mz, rrit, Mh', Ga is calculated. 



It is recommended to use the IM0MS=1 option. Other options are introduced for some tests. 

There are also some additional parameters which could be set by the user. For the 
TERAD branch one can set parameters GG , GZ , ZZ (=0 , 1) which switch on/off the contri- 
butions due to 7-exchange, 7Z-interference, Z-exchange, respectively, in cross sections in 
any variable in the TERAD branch. (Together with IN0RM=2 one must fix GG=GZ=ZZ=1.) 

There are two parameters W2TR and TTR used to subdivide the 2-dimensional inte- 
gration region in all DXY. . .'s into up to three parts. For W2TR=TTR=0, there is only one 
region. Choosing the parameter TTR (in GeV^) > 0, one splits the region into two parts, 
one of them with small Q\ where the structure functions should have a special behaviour. 
In addition, this region of small Q\ can be split into two subregions by a choice of the 
parameter W2TR (in GeV^) > Mp^^^^^. The latter was used for a proper treatment of the 
resonance region. Such subdivision of the kinematical region is used for MODEL=0. For 
more details about the low x region see [15]. 

Two more parameters are: ALAM - the lepton beam polarization (in the interval 
( — 1,+!)), and LEPCH - the lepton beam charge, —1 for electron, +1 for positron beams 
respectively. 



In SUBROUTINE SET one sets also some constants like n, proton and electron masses, 
electromagnetic coupling constant a, Fermi coupling G^, and some combinations of them. 

Then goes a part which initializes the calculation of electroweak parameters. It needs 
an additional setting of flags (they are stored in the NPAR-array) and constants for initial- 
ization of the EW library DIZET. Since this library was originally designed for LEP-I 
physics some flags and parameters have not too much meaning at HERA, therefore for 
some of them we simply report the recommended value avoiding further explanation. 
IHVP 

says how to calculate the hadronic vacuum polarization (HVP) contribution to the run- 
ning aQED- 

IHVP=1: HVP of Jegerlehner, [17], 
IHVP=2: HVP is made of some effective quark masses, 
IHVP=3: HVP of Burkhardt et al., [18], 
IHVP=3: recommended value. 
IAMT4 

says how to apply some leading higher order corrections, connected with the top-quark, 
to Z, W widths and EW form factors: 
IAMT4=0: they are neglected, 

IAMT4=1,2,3: different variants of their treatment, 
IAMT4=3: recommended value. 
IQCD 

says how to apply QCD 0{aas) corrections: 
IQCD=0: they are neglected, 

IQCD=3: they are included approximately (applicable only at LEP-I), 
IQCD=4: they are calculated completely, 

IQCD=4: recommended value (with a subsequent redefinition, see below). 
IMASS 

IMASS=0: recommended value. 
lALST 

IALST=1: recommended value. 
IQCD3 
IQCD3=1: recommended value. 

The Z-boson mass AMZ, top-quark mass AMT, Higgs boson mass AMH and VARQCD = 
asiM"^) (=0.12 recommended value) should be also set for the subsequent call of the EW 
library, CALL DIZET. For IM0MS=1, DIZET calculates the l^-boson mass AMW, but having 
in mind IMOMS 7^ 1 options, the W^-mass can also be set. In the ZPAR-array, DIZET returns 
many useful static EW quantities. We mention only two: Ar is stored in ZPAR(l) and 
w 



sin^ 6'w in ZPAR(3). 



After the DIZET call one could re-initialize, if needed, the EW parameters. To do that 
it is sufficient to redefine AMZ and AMW. There is then one more important redefinition. 
As the IQCD=4 option is very time consuming, we recommend to retain after the call of 
DIZET two statements: 
IQCD = 
NPAR(3) = IQCD 
which will result in applying 0{aas) corrections only in Ar but not in the EW form 



factors. For applications at HERA this is a very good approximation. The subsequent 
CALL CONSTQ initiahzes all fermion (lepton and quark) masses and charges. 

Then one sets some other constants like QE = 2J3, twice the weak isospin of the 
electron, GV and GA, vector and axial vector electron couplings to the Z-boson, and the 
pion threshold W2PIT = {Mproton + ^pion)^- 

SUBROUTINE SET ends with the calculation of the s-invariant at HERA and the cre- 
ation of a lattice in x and y over the kinematic domain in which DIS cross sections will 
be calculated. The version described creates the 6(a;)*16(|/)=96 points defined by 



x = .0001, .001, .01,.!, .5, .9 

and 

y = .01, .02, .05, .1, .15, .2, .3, .4, .5, .6, .7, .8, .9, .95, .98, .99. 



5 Further developments of the code 

The code TERAD91 was frozen in October '91 and used to produce numbers for various 
Figures and Tables of the Workshop Proceedings. Last bug fixes were made on February 
10th of 1992, for this reason it is recommended to use a version saved after that date. 
Nowadays only one new calculation is in progress: the DZRCG creates a DISEP branch 
for the NC reaction in mixed and hadronic variables. Later it is planned to add a DISEP 
branch for the CC reaction in hadronic variables. Further efforts in the future are not 
excluded if some new results will be requested by HERA experiments. 



The TERAD91 package consists of about 11150 lines of coding (including about 4000 
lines of weak library DIZET). 
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